library(limma)
library(affy)
library(geneplotter)

setwd("GLDilution")
pd <- read.phenoData("pdata.txt",as.is=TRUE)
sns <- paste("L","S_",pd$LiverAmt,"_",pd$SN19Amt,sep="")
fns <- paste(pd$filename,"gz",sep=".")

Data<-ReadAffy(filenames=fns[11:15],phenoData=pd[1:5,],description="miame.txt")
sampleNames(Data) <- LETTERS[1:5]

bitmap("../figure-01.png",res=300,pointsize=20)
mypar(1,1)
boxplot(Data)
dev.off()

bitmap("../figure-02.png",res=300,pointsize=20)
mypar(1,1)
hist(Data)
dev.off()

###Affy SpikeIn
library(SpikeIn)
data(SpikeIn95)
for(i in 1:2){
  if(i==1) Data <- SpikeIn95[,9:10] else Data <- normalize(Data)
  x <- log2(pm(Data[,1]))
  y <- log2(pm(Data[,2]))
  if(i==1){
    bitmap("../figure-03-0.png",res=300,pointsize=20)
    mypar(1,1)
    smoothScatter(x,y,xlab="Expression 1",ylab="Expression 2")
    dev.off()
  }
  A <- (x+y)/2
  M <- y-x
  Index <- which(abs(M)<1.2)
  A <- A[Index]
  M <- M[Index]
  bitmap(paste("../figure-03-",i,".png",sep=""),res=300,pointsize=20)
  mypar(1,1)
  smoothScatter(A,M)
  sIndex <- which(probeNames(Data)[Index]=="40322_at")
  points(A[sIndex],M[sIndex],col="red",pch=16,cex=.25)
  fit1 <- loess(M~A,subset=order(A)[seq(1,length(A),len=2000)],
                degree=1,span=1/3)
  if(i==1) simfit <- fit1 ###for next plot
  lines(sort(fit1$x),fit1$fitted[order(fit1$x)],col=2)             
  dev.off()
}

##Loess demo
A <- seq(min(simfit$x),max(simfit$x),len=2000)
M <- predict(simfit,A) + rnorm(length(A),0,0.15)
bitmap("../figure-04.png",res=300,pointsize=20)
mypar(1,1)
plot(A,M,cex=.25,col=3)
fit1 <- loess(M~A,degree=1,span=1/3)
lines(sort(fit1$x),fit1$fitted[order(fit1$x)],col=2)             
dev.off()
for(i in 1:5){
  bitmap(paste("../figure-04-",i,".png",sep=""),res=300,pointsize=20)
  mypar(1,1)
  plot(A,M,cex=.25,col=3)
  fit1 <- loess(M~A,degree=1,span=1/3)
  cutoff <- (8:12)[i]
  a <- A[which(A>cutoff-1 & A<cutoff+1)]
  m <- M[which(A>cutoff-1 & A<cutoff+1)]
  points(a,m,cex=.35,col=8,pch=16)
  lines(sort(fit1$x),fit1$fitted[order(fit1$x)],col=2,lwd=2)
  abline(lm(m~a),col=4,lwd=2,lty=2)
  dev.off()
} 

###beta7
setwd("../beta7")
targets <- readTargets("TargetBeta7.txt")

RG <- read.maimages(targets$FileName, source="genepix")
RG$printer <- getLayout(RG$genes)
types <- readSpotTypes()
RG$genes$Status <- controlStatus(types, RG)
MA <- MA.RG(RG,bc.method="none")

bitmap("../figure-05.png",res=300,pointsize=20,width=12)
mypar(1,1)
boxplot(split(MA$M[,2],MA$genes$Block),range=0,ylim=c(-2,2),xlab="Print-tip",ylab="M",col=rep(1:12,rep(4,12)))
dev.off()

bitmap("../figure-06.png",res=300,pointsize=20,width=12)
mypar(1,1)
imageplot(MA$M[,2], RG$printer, zlim=c(-3,3))
dev.off()

bitmap("../figure-07.png",res=300,pointsize=20,width=12)
mypar(1,1)
imageplot(log2(RG$G[,2]), RG$printer)
dev.off()

Order <- printorder(RG$printer)
bitmap("../figure-08.png",res=300,pointsize=20,width=12)
mypar(1,1)
boxplot(split(MA$M[,5],Order$plate),range=0,ylim=c(-2,2),col=rainbow(61))
dev.off()

x=Order$printorder
y=MA$M[,5]
bitmap("../figure-09.png",res=300,pointsize=20,width=9)
mypar(1,1)
smoothScatter(x[abs(y)<1.5],y[abs(y)<1.5],xlab="Print Order",ylab="M")
fit1 <- loess(y~x,degree=1,span=1/5,subset=sample(length(x),2000))
lines(sort(fit1$x),fit1$fitted[order(fit1$x)],col=2)             
abline(h=0,lty=2)
dev.off()

bitmap("../figure-10.png",res=300,pointsize=20)
mypar(1,1)
plotPrintTipLoess(RG)
dev.off()


x <- log2(100+2^seq(0,15,len=100))
y <- log2(50+2^seq(0,15,len=100))
A <- (x+y)/2
M <- y-x
bitmap("../figure-11.png",res=300,pointsize=20)
mypar(1,1)
plot(A,M,type="l",ylim=c(-1.2,1.2))
dev.off()

##f <- function(x) as.numeric(x$Flags > -75)
##RG <- read.maimages(targets$FileName, source="genepix", wt.fun=f)
##RG$printer <- getLayout(RG$genes)
##RG$printer
